function orientation_ch_norm = orientation_channel_sim(img, pars)

%% input parameters
img_size_x = pars.img_size_x;
img_size_y = pars.img_size_y;
Nc = pars.Nc;
Ns = pars.Ns;
delta = pars.delta;
c = pars.c;
pyr_filt_coeff = pars.pyr_filt_coeff;
resample_filt_coeff = pars.resample_filt_coeff;
levels = pars.levels;
out_scale = pars.outscale;
dec_factor = pars.dec_factor;
N_theta = pars.N_theta;
theta = pars.theta;

output_nr = img_size_y/(dec_factor^out_scale);  %size of output conspicuity map
output_nc = img_size_x/(dec_factor^out_scale);

cs_ch = zeros(Nc*Ns,2);
orientation_ch = zeros(output_nr,output_nc);

orientation_ch_acc = zeros(output_nr,output_nc);
% the center-surround levels for each channel
for i = 1:Nc
    clev = c(i);
    for j = 1:Ns
        slev = clev + delta(j);
        k = (i-1)*Ns+j;
        cs_ch(k,1) = clev;
        cs_ch(k,2) = slev;
    end
end

t_index = 1;
%% computation begins
% for t_index = 1:N_theta   %loop over theta
    
    pyr_select = 1; % 0 : gaussian pyramid, 1 : laplacian pyramid
    pyr = pyramid_sim(img, pyr_filt_coeff, levels, pyr_select);         %Laplacian Pyramid generation

     for pyrlev = 1:levels+1            %generate gabor pyramid
        pyr_gabor{pyrlev} = gabor_sim(pyr{pyrlev}, theta(t_index), pars);
     end
    
    for lpi = 1:Nc*Ns   

        clev = cs_ch(lpi,1); %center channel for csd
        slev = cs_ch(lpi,2); %surround channel for csd
        cs_factor = 2^(slev-clev);

        center = pyr_gabor{clev+1};   %+1 to for matlab array index
        surround = pyr_gabor{slev+1};

        csd_out = csd_sim(center, surround, cs_factor);     %Center-surround difference
     
     %norm_out = csd_sim(center, surround, cs_factor);

       norm_out = maxnorm_sim(csd_out);                    % Max normalizer

        resample_lev = abs(out_scale-clev);  % num of levels to output scale
        resample_out = resample_sim(norm_out, resample_filt_coeff, resample_lev); %Resample to output scale

        orientation_ch = orientation_ch + resample_out;    %Across-scale addition
    end

%     orientation_ch_theta = maxnorm_sim(orientation_ch);      %Final normalization for each theta
     orientation_ch_theta = orientation_ch;      
     orientation_ch_acc = orientation_ch_acc + orientation_ch_theta;  % Accumulate the orientation map for each theta
% end

orientation_ch_norm = maxnorm_sim(orientation_ch_acc);      %Final normalization to get the consicuity map
% orientation_ch_norm = orientation_ch_acc;     



